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ABSTRACT 


This document provides computational models for the Flight 
Experiment Demonstration System (FEDS) , developed at the 
Goddard Space Flight Center System Technology Laboratory 
(STL), Code 580. FEDS is a modification of the Automated 
Orbit Determination System developed at the STL during 1981 
and 1982. 
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SECTION 1 


INTRODUCTION 


This document provides computational models for the Flight 
Experiment Demonstration System (FEDS ) , developed at the 
Goddard Space Flight Center (GSFC) System Technology Labora- 
tory (STL ) , Code 580. FEDS is a modification of the Auto- 
mated Orbit Determination System (AODS ) , developed at the 
STL during 1981 and 1982. The purpose of FEDS is to demon- 
strate, in a simulated spacecraft environment, the feasibil- 
ity of using microprocessors to perform onboard orbit 
determination with limited ground support. 

During the demonstration, FEDS will execute on the STL's 
PDP-11/23 microcomputer, located at a remote site. FEDS 
will be connected to a transponder via a communications link 
and to the AODS Environment Simulator for Prototype Testing 
(ADEPT) , modified to support FEDS. The transponder will 
provide observation data, accumulated in real time from a 
Tracking and Data Relay Satellite System (TDRSS) signal, to 
FEDS for use in the orbit determination process. ADEPT, 
executing on the PDP-11/70, will provide all other external 
information required by FEDS. 

The models presented in this document were largely taken 
from the AODS System Description, Section 5 (Reference 1) . 

In most cases, only the algorithms are presented here; no 
attempt is made to derive models or show how they are inte- 
grated into FEDS. The reader should refer to the AODS Sys- 
tem Description, Section 3, for a logical description of the 
corresponding tasks. 

The coordinate systems, transformations, and time systems 
used throughout FEDS are described in Section 2. The orbit 
propagation algorithms, including integration and interpola- 
tion algorithms, the equations of motion, and the variational 
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equations, are described in Section 3. Algorithms for 
reduction of raw observation data and production of 
frequency control words are given in Section 4, and the 
estimation algorithms and observation models are presented 
in Sections 5 and 6, respectively. 
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SECTION 2 - COORDINATE SYSTEMS, TRANSFORMATIONS, AND 

TIME SYSTEMS 


This section contains a description of the coordinate sys- 
tems, transformations, and time systems used in FEDS. It 
also contains the algorithm for the Greenwich hour angle 
calculation, which is used in transformations throughout the 
system. 

2.1 COORDINATE SYSTEMS AND TRANSFORMATIONS 

In FEDS, the propagation of the satellite's state vector and 
state transition matrix is performed in geocentric rec- 
tangular coordinates referenced to a true Equator and 
equinox of epoch (TOE) coordinate frame. The satellite's 
acceleration vector and the partial derivatives of the ac- 
celeration vector are also expressed in this geocentric sys- 
tem. This system is obtained by freezing the true Equator 
and equinox of date (TOD) system at a specified epoch (ref- 
erence time) . The tracking measurements are also computed 
in the TOE system. However, the coordinates of the ground- 
fixed antenna are expressed in Earth-fixed coordinates and 
must be transformed to the TOE system. Each of these coor- 
dinate systems is defined in the following subsections along 
with the required transformations. 

2.1.1 TRUE EQUATOR AND EQUINOX OF EPOCH (INERTIAL) SYSTEM 
The geocentric TOE system is defined as follows; 

Origin Center of the Earth 

Reference plane Equatorial plane of the Earth, 

perpendicular to the Earth's spin 
axis at epoch 

Principal direction True vernal equinox of epoch 

The equinox is defined as the intersection of the planes of 
the Earth's Equator and the ecliptic. The Equator is de- 
fined as being normal to the Earth's instantaneous spin 
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axis. The ecliptic is the Earth-Sun orbital plane. The TOE 
system is equivalent to the TOD system associated with the 
specified epoch date. However, the TOD system changes with 
time because of the motion of the equinox, which is due, in 
turn, to the combined motions of the two planes (the Equator 
and the ecliptic) that define it. 

The motion of the Earth's spin axis or of the Earth's Equator 
is due to the gravitational attraction of the Sun and the 
Moon on the Earth's equatorial bulge. This motion consists 
of lunisolar precession and nutation. The motion of the 
ecliptic is due to the planets' gravitational pull on the 
Earth and consists of a slow rotation of the ecliptic. This 
motion is known as planetary precession. 

The rectangular Cartesian coordinates (see Figure 2-1) asso- 
ciated with the TOE coordinate system are defined with re- 
spect to the following axes: 

x-axis = principal direction 

y-axis = normal to the x and z axes to form a right-hand 
system 

z-axis = normal to the equatorial plane of epoch in the 
direction of the Earth's spin axis 

EARTH'S 
SPIN AXIS 



VERNAL EQUINOX n 

AT EPOCH « 


Figure 2-1. Geocentric Coordinate System 
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Quantities r, x, y, and z designate the position vector and 
Cartesian coordinates referred to the TOE frame. 

2.1.2 EARTH-FIXED SYSTEM 

The Earth-fixed coordinate system is defined as follows: 

Origin Center of the Earth 

Reference plane Equatorial plane of the Earth, 

perpendicular to the adopted 
polar geographic axis 

Principal direction Intersection of the Greenwich 

meridian with the reference plane 

The Earth's axis of figure (i.e., principal moment of in- 
ertia) is not coincident with the Earth's instantaneous spin 
axis. It moves with respect to the latter, causing the 
polar motion effect. Therefore, the motion of the spin axis 
pole is given with respect to the pole at some established 
epoch. The pole at the established epoch is referred to as 
the adopted geographic pole and corresponds to the Earth- 
fixed z-axis, The adopted geographic pole is the mean 

pole of 1903.0, which is consistent with that used by the 
International Polar Motion Service. 

The Greenwich meridian is the plane containing the adopted 
polar axis that passes through the former Royal Observatory 
at Greenwich, England. 

The rectangular Cartesian coordinates (see Figure 2-2) asso- 
ciated with the Earth-fixed coordinate system are defined 
with respect to the following axes: 

x^-axis = principal direction 

y^-axis = normal to the Xj^ and z^ axes to form a right- 
hand system 

zi 3 ,-axis = axis along the vector passing through the 
adopted geographic pole 
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ADOPTED 

GEOGRAPHIC 

POLE 



2.1.3 EARTH-FIXED-TO-TOE TRANSFORMATION 

The exact transformation that relates the TOE coordinates to 
the Earth-fixed coordinates accounts for three separate ef- 
fects. The first relates the true vernal equinox of epoch 
to the true vernal equinox of date (i.e., accounts for pre- 
cession and nutation effects) . The second effect relates 
the true equinox of date of the Greenwich meridian of the 
rotating Earth by means of the angle, the true right 

ascension of Greenwich (see Figure 2-3) . The third effect, 
called polar motion, accounts for the fact that the pole of 

the Earth-fixed axis, z. , does not coincide with the 

b 

Earth's instantaneous spin axis, the pole of the TOD geo- 
centric axes. 

In FEDS, the first and third effects are assumed to be 

zero. This is equivalent to assuming that the TOE z-axis is 

coincident with both the instantaneous spin axis and the 

Earth-fixed z, -axis. The transformation from the TOE to 
b 

the Earth-fixed coordinate system reduces to a rotation 
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about the TOE z-axis through the true right ascension of 

Greenwich, a^, yielding 
3 




( 2 - 1 ) 


where 


= 


cos a 


-sin a. 


sin a 0 

3 

cos ttg 0 


( 2 - 2 ) 


Computation of the true right ascension of Greenwich, a^, is 
discussed in Section 2.1.4. 


b GREENWICH 



(V 

CO 

*s» 

CM 


Figure 2-3. TOE-to-Ear th-Fixed Coordinate Rotation 


The Earth-fixed coordinates are transformed into TOE co- 
ordinates as follows: 



(2-3) 
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Differentiation yields the velocity transformation 


t,T 




where 



■-sin a 

g 

-cos Og 


0 


cos a 

g 

-sin Og 
0 



(2-4) 


(2-5) 


and where Og (the rate of change of ttg) is considered to be 
constant. 

2,1.4 COMPUTATION OP THE GREENWICH SIDEREAL TIME AND HOUR 
ANGLE 

The TOD right ascension of Greenwich, a , is measured 

y 

easterly from the true vernal equinox to Greenwich. A 
related quantity is the Greenwich hour angle (GHA) , also 
called the true Greenwich sidereal time, which is measured 
westerly in the plane of the Equator from Greenwich to the 
true vernal equinox. Thus, although their definitions 
differ, the right ascension of Greenwich, Og, and the 
Greenwich sidereal time and hour angle are equal in magni- 
tude. In FEDS software, the approximation is made that the 
TOE Greenwich sidereal time is equal to the TOD Greenwich 
sidereal time. 

The true Greenwich sidereal time, Ug, is obtained from the 
mean Greenwich sidereal time, by applying a correction, 

AH, due to nutation in longitude and obliquity as follows: 


a. 


“gm 


( 2 - 6 ) 
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The mean Greenwich sidereal time is calculated using the 
following equation; 


= 100?075542 + 0?98564735 Ad 
+ (2°9015 X lO"^^) Ad^ + 


(2-7) 


where is the angular velocity (degrees per second) of the 
Earth and is given by 


0?0041780742 

1.0 + (5.21 X 10‘^^) Ad 


( 2 - 8 ) 


where Ad = number of whole days elapsed from atomic time 
(A.l) January 0, 1950 

f (3 = fraction of days measured from midnight (seconds) 

If the input request time is expressed in terms of the A.l 
modified Julian date (i.e., days since 2430000), DJD, then 


DJO - 3282.5 = Ad + f , 

d 


(2-9) 


The nutation term, AH, is computed by evaluating a ninth- 
order nutation polynomial derived from the GTDS solar- 
lunar-planetary ephemeris (SLP) file (Reference 2, 

Section 3.6) 


10 

AH = ?! + ^ Pi(At)^"^ (2-10) 

i=2 


where P^ = coefficient of the ith term of the nutation 

polynomial that covers the request time, DJO 

At = DJO - (PDELHT + 10) (days) 
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PDELHT = uplink start of polynomial fit (modified 
Julian date in A.l days from start of SLP 
file) 

Since each nutation polynomial covers a 20-day span, two 
polynomials are used in FEDS to provide for continuity over 
the estimation data arc. Only the observation modeling rou- 
tines require the addition of the nutation term described 
earlier. 

Subsequent requests to compute at another time, t, may 

use a reference time and an associated a as follows: 

9ref 

“g = 

where a = reference GHA computed at time t^.^^ 

^ref 

(0 = rotation rate of the Earth 
2.2 TIME SYSTEMS 

The basic time system used in FEDS is the A.l system. Be- 
cause time advances at a constant rate in A.l time, simula- 
tion time can be compared directly with system clock time 
(in seconds from reference) for system control and schedul- 
ing purposes. All input time tags are in the Universal Time 
Coordinated (UTC) system and are converted to A.l time on 
input as follows: 


A.l 


"IN 


+ a 


il 


where 


t^^l = internal time (A.l) 
tjN = input time (UTC) 
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ail “ constantl adjustment (UTC to A.l) for time period 
surrounding tjfj 

After output, all time tags are converted from A.l to UTC as 
follows ; 


^OUT *^A.l " ^il 


(2-13) 


where t^.i = internal time (A.l) 
touT “ output time (UTC) 

ail “ constant^ adjustment (UTC to A.l) for time 
period surrounding tQjj.p 

One other time system is used in FEDS. The A.l time is con- 
verted to ephemeris time (ET) during the computation of the 
Sun and the Moon positions in the orbit propagator as 
follows ; 


’^ET ^A.l 86400 
2. 3 TIME REPRESENTATIONS IN FEDS 

Time is represented in three basic ways in FEDS. When the 
START command is received, a time message is requested from 
an external Parallel Grouped Binary Time Code 5 (PBS) gener- 
ator. This message is used to establish a simulation refer- 
ence time. All input time tags are converted from the input 
form of YYMMDDHHMMSS.SS (UTC) to seconds from reference 
(A.l). All time parameters are kept and measured in FEDS in 
terms of seconds from reference. System clock time is also 
measured in seconds from reference for controlling real-time 
processing . 


^The standard conversion polynomial is t^.i = ^li 

a 2 i At + a 3 i At^. However, since a2i and a 3 i are zero for 
the times FEDS will use, they were omitted to save memory. 
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The PBS time code used in FEDS comprises three groups, 
totaling 41 bits: truncated Julian's Day (TJD) , seconds of 

day, and milliseconds. TJD is defined as the number of 
whole days since midnight (UTC) of May 5, 1968, and is com- 
puted as follows: 

TJD = Julian Day Number - 244000.5 

The T JD-to-calendar date conversion is presented in 
Table 2-1. 

Modified Julian dates are also used throughout FEDS. A 
modified Julian date is defined as follows in FEDS: 

DJO = days since 2430000 in A.l time 

The reference time is also stored in modified Julian date 
form for quick computation of days since reference time. 
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Table 2-1. TJD-to-Calendar-Date Conversion 


TJD Date (Yeary Month, 


0 

680524* 

222 

690101 

587 

700101 

952 

710101 

1317 

720101* 

1683 

730101 

2048 

740101 

2413 

750101 

2778 

760101* 

3144 

770101 

3509 

780101 

3874 

790101 

4239 

800101* 

4605 

810101 

4970 

820101 

5335 

830101 

5700 

840101* 

6066 

850101 

6431 

860101 

6796 

870101 

7161 

880101* 

7527 

890101 

7892 

900101 

8257 

910101 

8622 

920101* 

8988 

930101 

9353 

940101 

9718 

950101 

0 

951010 

• 

• 

• 

• 


*Denotes Leap Year 
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SECTION 3 - ORBIT PROPAGATION MODELS 


Orbit propagation is performed in FEDS using a multistep 
integration process (with a single-step starter) . The 
propagator integrates one full fixed-length step at a time 
and uses a multistep interpolator to obtain the spacecraft 
state or state transition matrix at the request time. 

Because the FEDS orbit propagator is required to propagate 

high-altitude geosynchronous orbits (TDRS) and low-altitude 

drag-perturbed orbits (user spacecraft) , various forces are 

modeled in the equations of motion, including atmospheric 

drag, up to a 15-by-15 geopotential field, gravitational 

effects of the Earth, the Sun, and the Moon, and solar 

radiation pressure. However, due to timing and memory 

constraints, only the gravitational effect of the Earth and 

2 onal coefficients, X_ , X , and X_ , are used in the 

J2 J3 J4 

variational equations. 

The computational models for the multistep integrator and 
interpolator, the single-step integrator (starter) , the 
equations of motion, and the variational equations are given 
in the following subsections. 

3.1 RUNGE-KUTTA STARTER 

The multistep method avoids the multiple function evalua- 
tions at each integration step that are required by the 
Runge-Kutta single-step method, but it is not self- 
starting. Starting from an initial position and velocity, 
the Runge-Kutta method is used to build the required start- 
ing array for the Cowell equations of motion and variational 
equations. After the Runge-Kutta integrator has entered 
10 equally spaced accelerations in the backpoints table, the 
initial ordinate first and second sums are computed. The 
Runge-Kutta algorithm and the initial computation of the 
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first and second ordinate sums are presented in the follow- 
ing subsections. 

3.1.1 RUNGE-KUTTA INTEGRATION ALGORITHM 

The Runge-Kutta method is a single-step, self-starting nu- 
merical integration technique by which the value of the 
dependent variable at some future time, t^^ + x, where x is 
the starter integration step size (x = h/6, where h is the 
multistep integration step size) , can be calculated from a 
weighted summation formula and the value of the dependent 
variable at t^. The integrator algorithm is as follows. 

Given a first-order differential equation of the form 


il = 


X) 


and an initial value 


= x(t^) 

the dependent variable, x, at time t^^ + x is as follows: 

N 

x(tj^ + X) = x^ + X Wjfj (3-1) 

j=l 

where 

f^ = f(t^, x^) (3-2) 


3-2 


9681 



and 


f . 
D 


f 



j-1 
i=l 


(3-3) 


The derivatives f are evaluated at n subintervals from t^^ 
to t^ + X. The location of these subintervals is defined by 
coefficients Yj* The summation is weighted by and w j . 
The number of derivative evaluations needed is N, and the 
error made in the integration is of the order of x^^^, where 
p is the order of the method, which is not necessarily equal 
to N. 

Runge-Kutta methods have been developed by comparing two in- 
tegrators of different order in which the lower order method 
derivatives are a subset of those used for the higher 
order. An estimate of the truncation error made from using 
the lower order integration method can be obtained by dif- 
ferencing the two solutions of orders p and q (where 
q < p) to obtain 


N 

E(t^ + X) = X (Sjf j (3-4) 

j=l 

where N is the number of derivative evaluations for the 
solution of order p. Equation (3-4) can be used to compute 
the relative truncation error, at each integration 

step as follows; 


E(tj^ + X) 

lx ( tj^ + x) I 


(3-5) 
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Runge-Kutta integration formulas are available for which the 
integration constants have been chosen so as to minimize the 
truncation error of the higher order method and to use all 
free parameters in the equations of condition to optimize 
the integrator. This approach produces better integrators 
in that there are fewer function evaluations for the ac- 
curacy obtained than for Runge-Kutta integrators derived 
without such considerations. 

The integrator used in FEDS is a fixed-step integrator de- 
rived under these conditions. The integration coefficients 
(see Table 3-1 drawn from Reference 3) require five deriva- 
tive evaluations and produce a solution accurate to order p, 
where p is greater than four but less than five. The error 

coefficients, 6., are consistent with the difference of 

3 

solutions of orders p and q, where q equals three. For this 
reason, the algorithm is referred to as a Runge-Kutta 3(4+) 
method. 

3.1.2 INITIAL COMPUTATION OF THE ORDINATE FIRST AND SECOND 
SUMS 

The equations for computing the initial ordinate first and 
second sums for state integration are as follows: 


10 


^s = _n - V 6* X 

^n-1 h “^i+l ^n-i 


i+1 


(3-6) 


i=0 


10 


II, 


n-1 ^2 Lu i+2 n-i 


+1 


(3-7) 


i=0 


where n = 10 

Xj^ and Xj^ = position and velocity vectors at time tj^ 
integration 


3-4 


9681 



Table 3-1. Coefficients for Runge-Kutta 3(4+) Integrator 


cn/tcav 
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h = integration step size 

X = accelerations stored in the backpoints 
table 

= Adams -Moulton correction coefficients 
= Cowell correction coefficients 

These equations are derived by inverting the Adams-Moulton 
and Cowell correction formulas given in Section 3.2 (Equa- 
tions (3-15) and (3-17)). 

The following equations are used to update the sums to time 


n 


n-1 


+ X 


n 


(3-8) 





(3-9) 


The equations for computing the initial ordinate first and 
second sums for the partial derivative integration are as 
follows : 


where 


10 


n-1 




B* y 
•"i+l n+l-i 


i=0 




^n+l-i 


(3-10) 


(3-11) 


n = 10 

Y = partial derivatives of the velocity at t^ 

Y = partial derivatives of the position at 
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Y = partial derivatives of the accelerations 
stored in the backpoints table 

h, = the values defined previously 

These equations are derived by inverting the Adams -Moulton 
and Cowell correction formulas given in Section 3.2.1. 

The following equations are used to update the sums to time 


1 + Y (3-12) 

n n-i 

= ^’^Pn 1 + ^Pn (3-13) 

n n-i n 

3.2 MULTISTEP INTEGRATION 

Multistep integration methods minimize the number of deriva- 
tive evaluations required to produce a given accuracy at the 
end of the requested interval of integration. Since, in 
general, the major cost in computing an orbit is the evalua- 
tion of the complex force function (see Section 3.4), multi- 
step algorithms appear to be most efficient. FEDS uses a 
12th-order, fixed-step, multistep integrator. A predict/ 
evaluate der ivative/correct/pseudoevaluate (PECE*) algorithm 
is used to integrate the equations of motion, and a 
corrector-only algorithm is used to integrate the varia- 
tional equations. The summed ordinate predictor/corrector 
formulas are given in Section 3.2.1, and the algorithms for 
integrating the equations of motion and the variational 
equations are described in Sections 3.2.2 and 3.2.3, respec- 
tively. Finally, the algorithm for computing the state 
transition matrix for mapping the partial derivatives is 
given in Section 3.2.4. 
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3.2.1 PREDICTION AND CORRECTION FORMULAS 


In FEDS, the Adams-Bashf orth predictor and the Adams -Moulton 
corrector methods are used to integrate first-order equa- 
tions of motion and the velocity vector. The Stormer pre- 
dictor and the Cowell corrector methods are used for 
second-order equations of motion of the position vector. 

The Adams-Bashforth prediction formula used to predict 
velocity is as follows: 


X 


n+1 


= h 



(3-14) 


where n is the number of backpoints available and k+1 is the 
order of the integrator. The Adams-Bashforth llth-order 
coefficients, are given in Table 3-2, 

The Adams-Moulton correction formula used to correct veloc- 
ity is as follows: 


k 

i=0 

where n is the number of backpoints in the table and k+1 is 
the order of the integrator. The Adams-Moulton correction 
coefficients, B*, i.e., the llth-order coefficients for equa- 
tions of motion and the 7th-order coefficients for varia- 
tional equations, are given in Tables 3-3 and 3-4, 
respectively. 
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Table 3-2. Adams-Bashforth llth-Order Predictor Coefficients 


= 3.451988400456282 
^2 = -13.81683726526174 
^3 = 35.43365334896585 
^4 = -60.83539672518839 
^5 = 72.18373924011945 
^6 = -59.70768131463444 
^7 = 33.93953914141414 
^8 = -12.67749704385121 
^9 = 2.808681814424002 
^lo = -0.2801895964439367 
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Table 3-3. Adams-Moulton llth-Order Corrector Coefficients 


n 

= 

0.2801895964439367 

6*2 

= 

0.6500924360169152 

6*3 

s 

-1.208305425284592 

65 

= 

1.810901775693442 

6*5 

= 

-1.99558147196168 

6*6 


1.575960936247395 

6*7 

= 

-0.8678660614077281 

6*8 

= 

0.3167875681417348 

6*9 

= 

-0.06896520387405804 

* 

10 

= 

0.006785849984634707 
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Table 3-4. Adams -Moulton 7th-Order Corrector Coefficients 


6| = 0.3155919312169312 
6* = 0.3921792328042328 
6^ = -0.3760251322751323 
8* = 0.2440806878306878 
8| = -0.09009589947089947 
3| = 0.01426917989417989 
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The Stormer prediction formula used to predict position is 
as follows; 


X 


n+1 



^n-1 


(3-16) 


where n is the number of backpoints in the table and k+2 is 
the order of the integrator. The 12th-order Stormer pre- 
diction coefficients, are given in Table 3-5. 

The Cowell correction formula used to correct position is as 
follows ; 

k 

" L “5+2 "n+l-i 

i=0 

where n is the number of backpoints in the table and k+2 is 
the order of the integrator. The Cowell correction 
coefficients, a^, i.e., 12th-order coefficients for equa- 
tions of motion and 8th-order coefficients for variational 
equations, are given in Tables 3-6 and 3-7, respectively. 

3.2.2 PECE* ALGORITHM FOR EQUATIONS OF MOTION 

The concept of pseudoevaluation is introduced as a device to 
help stabilize the numerical integration at little or no 
cost in computation. It is recognized that in a predictor- 
corrector scheme, the numerical stability region is propor- 
tional to the number of derivative evaluations within a 
given step (Reference 4) and that for systems of the form, 

X = f(X) + eg(X), where e is a small parameter, the 
stability region is mainly influenced by the f(X) term. 

The idea, then, is to introduce into a predictor-corrector 
algorithm, which is designed to solve the system above, a 
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Table 3-5. Stormer 



12th-0rder Predictor Coefficients 

0.7093330001402910 

-2.949775927679574 

7.565365635521886 

-12.95853327421036 

15.34411576078243 

-12.67150744799182 

7.193720363355780 

-2.684381939434023 

0.594237726972102 

-0.05924056412337662 
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Table 3-6. Cowell 12th-Order Corrector Coefficients 

a* = 0.05924056412337662 

a* = 0.1169273589065256 

ot* = -0.2839505421276255 
4 

a* = 0.4564979407166907 

a* = -0.5180148083012666 
6 

a* = 0.4154936016915184 
a* = -0.2309889820827321 
a* = 0.08485266855058522 
a*o = -0.01855655388207472 
a* = 0.001832085738335738 


9681 
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Cowell 8th-Order Corrector Coefficients 


Table 3-7. 


a* = 0.06820436507936508 

a* = 0.05115740740740741 

a* = -0.07000661375661376 
4 

ct* = 0.05019841269841270 
a* = -0.01936177248677249 
a* = 0.003141534391534392 


9681 


3-15 



pseudoevaluation, i.e., a partial evaluation of X, where 
f(X) is recomputed using the latest corrected value of X, 
and g(X) is reused based on a previous value of X. For ex- 
ample, assume that the equations to be integrated have the 
form 


*^ = ^ + P(t, r, t) 

r-^ 


(3-18) 


where the first term represents the primary attracting body 

acting on the satellite. Assuming that the accelerations 

and sums are known, then the iterative algorithm to advance 

to time t ,, is 
n+1 

1. Predict . Using Equations (3-14) and (3-16), predict 
values (denoted by superscript p) 


) = 




(p) 



n+1' 

^n+1' 

^n+lj 

(P) 

a(P) 


n+1' 

^n+1' 

^n+lj 


(3-19) 


(3-20) 


2. Evaluate. Using Equation (3-18) , evaluate 


-yr 




(p) 

n+1 


^(P)3 

^n+1 


•n+1' 


‘(P) 

n+1' 


-(P)\ 
n+1/ 


(3-21) 


3. Correct . Using Equations (3-15) and (3-17), obtain the 

— ^( c ) 

improved values (denoted by superscript c) , and 

i^(c) 

^n+1* 
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4. Pseudoevaluate. 


Compute the acceleration 


-yr 




■(c) 

n+1 


(c)3 


+ P 


n + 1 




-^(P) 

+1' ’^n+l' 


^(P)\ 

^n+1/ 


(3-22) 


where the P term is obtained from step 2. 
5. Update Sums . Compute the updated sums 


I 


n+1 



+ r 


^ ^n+1 ^ 


(3-23) 





(3-24) 


The computational cycle of steps 1 through 5 may then be 
repeated with n = n + 1. 

3.2.3 CORRECTOR-ONLY ALGORITHM FOR VARIATIONAL EQUATIONS 

In the Cowell formulation, the position and velocity partial 
derivatives of the satellite motion with respect to any pa- 
rameter appearing in the acceleration model in Equa- 
tion (3-18) or state (dynamic parameters) may be obtained by 
the numerical integration of a system of equations of the 
form 


Y = A(t) Y + B(t) Y + C(t) 


from initial conditions at time tg given by 


3r(t ) . 3r(t ) 

Y(t^) = — Y(to) = — ^ 


3p 


3p 


(3-25) 


(3-26) 
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where 


A(t) 


B(t) 


3r(t) 


3x3 


3r^(t) 

3 ^ 


3x3 


(3-27) 


(3-28) 


C(t) 


3r (t) 
a^p 


3 X Jl matrix of acceleration 
partial derivatives 


(3-29) 


Y(t) 


3r^(t) 

3P 


3x5, matrix of position partial 
derivatives 


(3-30) 


Y(t) 


3"^(t) 

3"p 


3x2, matrix of velocity 
partial derivatives 


(3-31) 


Vector p contains the parameters in the acceleration model 
to be estimated. 

The components of matrices A, B, and c are developed in Sec- 
tion 3.5. 

Optionally, the components of correspond to the space- 
craft's position and velocity at epoch and are expressed in 
true of date Cartesian coordinates. Since the first six 
elements of ^ are the state vectors, the first six columns 
of C are zero. The model parameter, drag, enters into 
P(t, r“, iT) of Equation (3-18) linearly, so that the computa- 
tion of C(t) is simplified by retaining many of the quanti- 
ties used in the computation of ^(t) . 

The integration of system Equation (3-25) is performed by 
utilizing the corrector-only formula as follows. Assume 
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that the satellite position and velocityr ^(^n+1^ ^ '^^n+1^ ' 

the matrices, step i, where i = 0, 1, 2, ..., k, 

and the summation matrices, and (3-by-il, are known; 

the algorithm to advance Y to time is 

1. Compute matrices and C(t^^j^), which 

depend only on and ‘^ri+l* 

2. Compute the 6-by-6 matrix [I - H]~^ 


‘'^“0*n+l 




H = 


^^O^n+1 '^^o\+l 


(3-32) 


where and 6* are the corrector coefficients of Equa- 
tions (3-15) and (3-17), and h is the step size. 

3. Form the 3-by-i. matrices, and as follows: 


X 


n 



^n+l-i *^0^n+l 


(3-33) 


V 


n 



n 


'n+l-i 




(3-34) 


4. Compute the required position and velocity partial de- 
rivatives, and by the matrix equation 


■^nn' 

, -1 

— 

^n 

• 

y , 

= 'I - »"6x6 

V 

n+1 

6xJl 

n 


(3-35) 


L J 6xil 
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5. Update acceleration and sums by 


= A Y + B Y + C (3-36) 

n+1 n+1 n+1 n+1 n+1 n+1 


I I •• 

P = P + Y 
n+1 n n+1 


(3-37) 


Hr, . I,, 

P , = P + P 1 1 
n+1 n n+1 


(3-38) 


and complete the cycle. After computing ^^id ^^^+ 2 ' 

steps 1 through 5 may be repeated with n = n + 1. 

3.2.4 STATE TRANSITION MATRIX COMPUTATION 

In FEDS, the orbit propagator computes (if requested) the 
state transition matrix that is used in the observation 
model to map the observation partial derivatives with re- 
spect to the solve-for parameters from the observation time 
to the epoch time. Since the orbit propagator is always 
restarted at epoch each time the epoch is changed, the state 
transition matrix is as follows: 





(3-39) 


where t = request time 
n 

t^ = epoch 

Y = partial derivatives of the position at time t^ 

with respect to the solve-for parameters (see 
Equation (3-30) ) 

Y = partial derivatives of the velocity at time tj^ 

with respect to the solve-for parameters (see 
Equation (3-31) ) 
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3.3 MULTISTEP INTERPOLATION 


The multistep interpolator uses the first and second ordi- 
nate sums and the 10 backpoints (accelerations) computed 
during multistep integration to compute the spacecraft posi- 
tion and velocity and the associated partial derivatives (if 
desired) at a request time. 

The multistep interpolation algorithm is as follows: 

1. Compute constants y]^(0) y^( 0) to be used in later 

calculations as follows: 


Yj(0) = 1 


y!(0) = 


i-1 

-r 

j = 0 


r-r-3-n Yj(0) 


(3-40) 


'■? (0) = ^ (0) Yj.jCO) 

j = 0 


(3-41) 


where 1=1, 2 , ..., 10. 
t - t 


2 . Le t s = 


n . 


in the following steps 


where t^^ = time associated with most recent entry in the 
backpoints table, Xn 

t = request time 

h = step size 

Compute Yj^(s) f where i = 0, 1, ..., 12: 


Yo(s) = 1 
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(i — If 2f , m m f 


12 ) 


(3-42) 


Yi(s) = 


s + i - 1 




4. Compute Y^(s), where i = 0, 1, 11: 


Y^(s) 


1 


Y-(s) = Y ■ (0) Yi_j (s) (i = 1. 2f ...f 11) 
j=0 

5. Compute 6|(s), where i = 0, 1, ...f 10: 


<S-(s) 


11 




Yi(s) 


(T) 


6^(s) = (-1)^ 

11 


* Z) r 

j=i+2 ' 

?) 

(i = 

1/ 2f •••f 9) 


is computed as 

m! 


(m-i) I i ! 



6'g(s) = Yl^(s) 


6. Compute the velocity as follows: 


x(t) = x(t^ + sh) 
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where = first ordinate sums 
n 

X = accelerations in the backpoints table 

n = number of accelerations in the backpoints 
table 

and s, h, and 6|(s) have been defined previously. 


7 . 


Compute Y^(s)f where i = 0, 1, 12; 


y;(s) = 1 


Y'!(s) = ^ Y- (0) 
j=0 


(i = 1, 2, 12) (3-46) 


8. Compute 6V(s), where i = 0, 1, 10: 


12 

<5;{s) = ^ y'^is) 
j=2 


6': (s) 


= (- 1 )^ 


12 


tV^i(s) + 




j=i+3 

(i = 1, 2, 9) 


Yj (s) 


(3-47) 




9. Compute the position as follows; 


X (t) = X (t^ + sh) = h' 


10 


- 1 ) 


i=0 


(3-48) 
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where and = first and second ordinate sums 

n n 

X = accelerations in the backpoints table 

n = number of entries in the backpoints 
table 

and s, h, and 6V(s) have been defined previously. 

The formulas for interpolating the partial derivatives are 
obtained by substituting 3x^_j^/3p for for and 

for in Equations (3-45) and (3-48) . 

3.4 EQUATIONS OF MOTION FOR STATE VARIABLES 

The orbital equations of motion, expressed in Cartesian co- 
ordinates, are 



(3-49) 


where r = satellite position vector in the TOE coordinate 
frame 

^ = total acceleration vector in the TOE coordinate 
frame 

In FEDS, this set of three second-order differential equa- 
tions is transformed to an equivalent set of six first-order 
differential equations; 



d^ 

dt 


(3-50) 


where r is the satellite velocity expressed in the TOE co- 
ordinate frame. 
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For the total FEDS state vector 


X = 


r 

r 

C, 


( 3 - 51 ) 


where Cj^ = atmospheric drag constant 


b = clock bias coefficients (b^^, b2f and b^) 


the equations of motion can be expressed as 


X = F(X, t) 


( 3 - 52 ) 


where 


F(X, t) = 


The total acceleration of the satellite includes the follow- 
ing components: 

• Gravitational acceleration of the satellite due to 
the Earth's mass C?_) and the solar and lunar 

t* 

masses (^g and respectively) 

• Gravitational acceleration of the satellite due to 
the nonsphericity of the Earth's gravitational po- 
tential (^5) 

• Satellite acceleration due to atmospheric drag 
forces (aj^) 
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• Satellite acceleration due to solar radiation pres- 
sure (a'gj^) 

The total acceleration, a, is expressed in terms of these 
components as 


a = a^ 


+ ag + 3 m + a 


*NS 


^SR 


(3-53) 


All or any subset of these effects can be included in the 
acceleration vector, which is used in constructing the equa- 
tions of motion. These accelerations are discussed in the 
following subsections. 

3.4.1 EARTH, SOLAR, AND LUNAR POINT MASS ACCELERATIONS 

To first order, the gravitational attraction of a body of 
mass m can be approximated as that arising from a dimension- 
less particle of mass m located at the center of mass of the 
body. Three point mass terms are included in the accelera- 
tion model of FEDS; the Earth's central body term and the 
perturbing accelerations from the Sun and the Moon. 

The acceleration experienced by a satellite attracted by n 
point masses, expressed in an inertial coordinate system, is 


df? _ 

dt2 = ■ rL 


""kp 


(3-54) 


k=l '■kp 


where r' = vector from the center of an inertial coordinate 
system to the satellite 

U)^ = product of the universal gravitational constant 
and the mass of the kth point mass 

= vector from the kth point mass to the satellite 

rkp = magnitude of the vector, 

In FEDS, the motion of the satellite is referenced to the 
Earth's position; i.e., an Earth-centered TOE coordinate 
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system is used. The desired form for the acceleration is 
obtained by subtracting the acceleration acting on the ref- 
erence body 


d r. 


dt' 


n 




k=l 


3 - ^k 


from each side of Equation (3-54) to obtain 


(3-55) 


d 2 . * 

(r - rj.) 

dt2 


n 




n 


=1 '•‘p 


kp 


-E 

k=l 

k?^E 


.3 


(3-56) 


For the case in which the reference central body is the 
Earth and the perturbing point masses are the Sun and the 
Moon, Equation (3-56) becomes 


where 
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The quantities Pg, and p^^ are the gravitational con- 

stants of the Earth, the Sun, and the Moon, respectively, 

and "rl and rt, are the position vectors of the Sun and the 
s M 

Moon, referenced to the Earth. 

The positions of the Sun and the Moon are determined by 
evaluating series expansions in the latitude and longitude 
of the Sun and the Moon, where the series are truncated from 
those of Brown's lunar theory by Woodard (Reference 5, 
pages 52 through 62). The Sun's position components are 
(Xg, yg, 2 g) , and the Moon's position components are (x^^, 

^M^ • 

For the Sun, the components are related to the orbital ele- 
ments by 


— = cos X, 


— = cos i sin X, 


(3-61) 


— = sin e sin X„ 

■'s ® 


For the Moon, the components are related to the orbital ele- 
ments by 




— - = cos cos X„ 
r.. M M 

M 


•M 


(3-62) 


= cos 9., sin X„ cos e - sin 0.. sin e 
r,, M M M 

M 


M . _ . _ 

— = cos 9„ sin X„ sin e + sin 0,, cos e 
r,, M M M 

M 
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where e = mean obliquity of the ecliptic 
Xg = solar longitude 
Xj^ = lunar ecliptic longitude 
= lunar ecliptic latitude. 

The distances r„ and r_ are found from evaluation of the 

M S 

terms dg/rg (the ratios of the mean to the true 

geocentric distances) as follows: 



1 + K 


M 



1 + Kg 


(3-63) 


where d^/I = mean lunar geocentric distance 

(384399.06 kilometers) 

dg = mean solar geocentric distance 
(149497971.0 kilometers) 

i<Mr " series presented in this subsection 

- 4 

The time AD (the elapsed time since 1900.00 in units of 10 
days) in the series expansion is given by 


AD = 3.6525 T 

c 


(3-64) 


where T^, the number of Julian centuries since 1900.00, is 
c 


JD - 2415020.0 
c ~ 36525.0 


(3-65) 


where JD is the full Julian date. 
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The mean obliquity of the ecliptic, e, is as defined as 
follows : 


e = 23?4522940 - 0?0130125 T 

c 

The solar longitude, Xg, is 

*S = ^ sin(lg) 

where 

Gg = 281?220833 + 4?70684 AD 
Hg = 358?475845 + 98560?0267 AD 
eg = 0.01675104 - 1.1444 x lO"^ AD 

In Equation (3-67) , Gg is the longitude of perigee of the 
solar orbit; Ag is the solar mean anomaly; and eg is the 
eccentricity of the solar orbit. 

The series k is given by 

O 

Kg = 6g cos(ilg) (3-69) 

The lunar ecliptic longitude, Xj^, is determined by subtract- 
ing a series in (Dj^, S,g) (presented in Table 3-8) 

from (the mean longitude of the Moon, measured in the 
ecliptic plane from the mean equinox of date to the mean 
ascending node of the lunar orbit, and then along the orbit) : 

^ ■ series 
3-30 
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Table 3-8. Series for Xn^ 


Coefficient 

(radians) 


"m 

Argument 

F 

Multiple of 

Dm 

M 


-0.000607 

s in 

0 

0 

1 

0 

0.11490 

sin 

0 

0 

2 

0 

-0.000267 

sin 

0 

2 

-2 

0 

-0.001996 

sin 

0 

2 

0 

0 

-0.000801 

sin 

1 

2 

-2 

0 

-0.003238 

sin 

1 

0 

0 

0 

-0.000118 

s in 

1 

0 

2 

0 

0.000138 

sin 

1 

0 

-2 

-1 

0.000716 

sin 

1 

0 

0 

-1 

0.000192 

sin 

1 

-2 

0 

0 

-0.000186 

sin 

1 

0 

-4 

0 

-0.022236 

sin 

1 

0 

-2 

0 

0.10976 

sin 

1 

0 

0 

0 

0.000931 

sin 

1 

0 

2 

0 

-0.000219 

sin 

1 

2 

0 . 

0 

-0.000999 

sin 

1 

0 

-2 

1 

-0.000532 

sin 

1 

0 

0 

1 

-0.000149 

sin 

2 

0 

-4 

0 

-0.001026 

sin 

2 

0 

-2 

0 

0.003728 

sin 

2 

0 

0 

0 

0.000175 

sin 

3 

0 

0 

0 
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The lunar ecliptic latitude, 9^^, is directly determined from 
a series in (Dj^, Fj^, J,g) and is presented in Table 3-9. 

The series is determined from the series in (D.., F.., 

M M M 

^M' presented in Table 3-10. In these series, 

= 270°26'o 2!'99 + (480960?0 52'59”31)T^ - 4"o8 T^ 

= 296°06'l6"59 + (477000°0 + 198°5o' ss” 79) T^ 

" 2 
+ 33.09 TZ 
c 

Ag = 358°28 '33 !'oO + (35640?0 + 359°02 ' 59 ”l0) T^ (3-70) 

= 11°15'03" 20 + (483120°0 + 82°0l' 3o" 54) T^ 

■I 2 

-11.56 T^ 

= 350°44 'l4l'95 + (444960?0 + 307° 06'5l"l8)T^ 

" 2 

- 5.17 T^ 

where Fj^ = argument of latitude of the Moon 

Djyj = mean elongation of the Moon from the Sun 

As noted previously, the series are summarized in Tables 3-8 
through 3-10. These tables give the coefficients of each 
trigonometric term and the angles that appear in that term. 
The trigonometric terms contain only integer multiples of 
the four angles. The first term in 6,^ is 

0., = 0.089503 sin(F„) radians 
M M 
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Table 3-9. Series of 0^ 


Coefficient 
( radians) 



Argument 

F 

M 

Multiple of 


0.089503 

s in 

0 

1 

0 

0 

0.000569 

sin 

0 

1 

2 

0 

-0.003023 

sin 

0 

1 

-2 

0 

-0.000144 

sin 

0 

1 

-2 

1 

0.004897 

s in 

1 

1 

0 

0 

-0.000807 

sin 

1 

1 

-2 

0 

0.004847 

s in 

1 

-1 

0 

0 

-0.000967 

sin 

1 

-1 

-2 

0 

0.00301 

sin 

2 

1 

0 

0 

0.000154 

sin 

2 

-1 

0 

0 

0.000161 

s in 

1 

-1 

2 

0 
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Table 3-10. Series for <{4 


Coefficient 

(radians) 



Argument 

Multiple of 


0.0082488 

cos 

0 

0 

2 

0 

0.0005604 

cos 

0 

0 

-2 

1 

0.0003369 

cos 

1 

0 

0 

-1 

-0.0002086 

cos 

1 

-2 

0 

0 

0.0100247 

cos 

1 

0 

-2 

0 

0.0545008 

cos 

1 

0 

0 

0 

0.0009017 

cos 

1 

0 

2 

0 

0.0004219 

cos 

1 

0 

-2 

1 

-0.0002773 

cos 

1 

0 

0 

1 

0.0029700 

cos 

2 

0 

0 

0 

0.0001817 

cos 

3 

0 

0 

0 
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The formulas for the lunar position include the leading 
terms of Brown's lunar theory as summarized by Woodard (Ref- 
erence 5) . All perturbation terms with amplitudes greater 
than 50 kilometers are included (22 terms in the ecliptic 
longitude, 11 terms in the ecliptic latitude, and 11 terms 
in the distance) ; this achieves an overall positional ac- 
curacy of 1 arc-minute (0.005 radian or 200 kilometers) . 

3.4.2 NONSPHERICAL GRAVITATIONAL ACCELERATION 

The inertial acceleration vector resulting from nonspherical 
gravitational effects is given by the gradient of the non- 
spherical terms in the geopotential function, 

^NS " ’^NS (3-71 


Expressed in terms of spherical polar coordinates, 

r = magnitude of the vector from the Earth's center 
of mass to the satellite 

(t) = geocentric latitude 

X = geocentric longitude (measured east from the 
prime meridian) 


The nonspherical 

geopotential. 

'I^NS" is given by 


00 

,,, > n 

't^NS ^ ^ ^ 



P°(sin (j>) 


n=2 




CO 

n 



^7^ Z 

z 

/R 

\-r) = 


n=2 

m=l 



* sin mX + c|!^ cos mxj 

where ]jg = gravitational constant of the Earth 

Rg = Earth's equatorial radius 
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P^(sin (t>) = associated Legendre function 

^n' *“n “ harmonic coefficients (zonal harmonics for 
m = 0 , sectorial harmonics for m = n, and 
tesseral harmonics for n > m 7 ^ 0 ) 

(NOTE: J_ = -C®, where J are the zonal 

n n 

coefficients 

The first and second terms are the nonspherical potential 
due to the sum of zonal and tesseral harmonics, respec- 
tively. The term n = 1 is not present, since the origin of 
the coordinate system is placed at the center of mass of the 
Earth. 

Expressing the gradient in Earth-fixed coordinates, 

= (X|^, yj^, Zj^) (see Section 2.1.2), the form for the in- 
ertial acceleration vector is obtained as follows; 



where x. / f ^nd z. are the components of the inertial 

°NS °NS °NS 

acceleration of the spacecraft, expressed in the Earth-fixed 
coordinate system and not the acceleration with respect to 
the Earth-fixed coordinate system. Thus, it is necessary to 
transform these components into the TOE coordinate system in 
which the orbital equations of motion are expressed. 
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This transformation, which is discussed in Section 2.1.2, is 
given by 



T 


NS 


(3-74) 


where the rotation matrix is defined in Equation (2-2) , 
and it is assumed that the geographic pole axis, z^, is 
aligned with the spin axis, z, of the TOE coordinate system. 
This rotation is equivalent to replacing (rj^, Xj^, z^) 

in Equation (3-73) by (r, x, y, z) , the TOE components, and 
calculating the longitude and latitude as follows: 


X = a - a. 


(fi = sin 


-1 


(!) 


(3-75) 


where a = right ascension of the spacecraft [a = tan”^ (y/x) ] 
a = right ascension of Greenwich 

g ^ 

Equation (3-73) can be written as 
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The partial derivatives of the nonspherical portion of the 
Earth's potential with respect to r, cj), and X are given by 
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(3-77) 


The Legendre functions and the terms cos mX, sin mX, and 
m tan t are computed via recursion formulas, as follows: 

P°(sin 4>) 

P^(sin (|)) 


P^(sin 4>) = (2n - 1) cos (j) (sin <()) (m ^ 0; m = n) 


(2n - 1) sin ^ P^_^^(sin (J>) 


- (n - 1) P°_ 2 (sin 4)) 

.m 


/n 


.m-1 


P^_ 2 (sin (])) + (2n - 1) cos P^_^(sin <t)) 
(m ^ 0; m < n) 


(3-78) 
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where 


P^Csin (p) = 1; P^(sin (|>) = sin <j); pj‘(sin <)>) = cos <J> 

sin mX = 2 cos X sin(m - 1) X - sin(m - 2) X 

cos mX = 2 cos X cos (m - 1) X - cos (m - 2) X 

m tan (j) = [ (m - 1) tan (j>) ] + tan cj) 

The partial derivatives of r, <j>/ and X, with respect to x, 
y, and z, are computed from the expressions 


3r 



dT r 


9 (|) 

dT 




+ y^) 


= 1 

^ 2 2 
aT (x^ + y^) 




(3-79) 


where the row vectors (1, 0, 0) r (0, 1, 0) , and (0, Of 1) f 
respectively, are given by 


3x ^ dz 

9 9 

dT dT a"? 
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Substituting Equations (3-79) into Equation (3-73) yields 
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(3-80) 


~T—2 

X + y 


''^NS 

3X 


‘NS 
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34 » 


NS 


3r 
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34» 


NS 


34> 


3.4.3 ATMOSPHERIC DRAG ACCELERATION 

Atmospheric drag acceleration is modeled as a drag force in 
the direction of the relative wind vector acting on a satel- 
lite of constant surface area. The velocity of the satel- 
lite relative to the atmosphere is computed in the inertial 
coordinate system by subtracting the motion of the atmos- 
phere, assumed to rotate with the Earth, from that of the 
satellite ; 



r - 0 ) X r 


(3-81) 


The Earth's rotation vector, w, is directed along the 
Earth's instantaneous spin axis with a magnitude equal to 
the rotation rate of the Earth and components (w^, 0 ) 2 , • 
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In the TOE coordinate system, ignoring the effects of pre- 
cession and nutation, the z-axis is aligned with the north 
polar spin axis such that u)^ and are equal to zero. In 
addition, o) is assumed to be constant. Therefore, in the 
TOE system. Equation (3-81) reduces to 



X + cjjy 

y - tox 
z 


(3-82) 


For the case of a spherical satellite, the atmospheric drag 
acceleration is computed as 


= - I (^) Oo'h) '%el' <3-® 

where Co = aerodynamic force coefficient, which- is an ad 
justable parameter 

A = surface area of the satellite 

m = mass of the satellite 

PQ(h) = altitude density function computed from the 
atmospheric drag model 

Nominally, for a spherical satellite, the aerodynamic force 
coefficient, is equal to 1.0. In order to absorb an 

error in any of the above terms, is an adjustable 
parameter . 

In FEDS, the altitude density function, Pp(h), is modeled 
using a Harris-Pr iester atmospheric model. Harris and 
Priester determined the physical properties of the upper 
atmosphere theoretically by solving the heat conduction 
equation under quasi-hydrostatic conditions (see Refer- 
ences 6 through 8) . Approximations for fluxes from the 
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extreme ultraviolet and corpuscular heat sources were in- 
cluded, but the model averages the semiannual and seasonal- 
latitudinal variations and does not attempt to account for 
the extreme ultraviolet 27-day effect. 

The atmospheric model presented here is a modification of 
the Harr is-Priester concept. The modification attempts to 
account for the diurnal bulge by including a cosine varia- 
tion between a maximum density profile at the apex of the 
diurnal bulge (which is located approximately 30 degrees 
east of the subsolar point) and a minimum density profile at 
the antapex of the diurnal bulge. Discrete values of the 
maximum- and minimum-density altitude profiles, correspond- 
ing to mean solar activity, are stored in tabular form as 

p (h.) and p . (h.), respectively. Different maximum and 
max 1 min i 

minimum profiles are available for different levels of solar 

activity (Reference 2). Exponential interpolation is used 

between entries; i.e., the minimum and maximum densities, 

p . and p , are given by 
min max' ^ ^ 


Ih ■ 

\ min / 

(3-84) 

Wh) = Omax*”!* ("B— ) 

\ max / 
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where (h^ £ ^ £ ^i+1^ respective scale heights, 

and are given by 

min max ^ 


H. 


"i - "i+i 


min 


in [p . {h . +1) /p . (h . ) ] 
min 1 min i' 


(3-85) 


H 


h . - h , . . 
1 1+1 


max 


in Ip (h . +1) /p (h . ) ] 
*^max ' 1 ' ^max i' 


A good approximation (neglecting polar motion) for the 
satellite height, h, is given by 


h = r - 


(3-86) 


where is the mean radius of the Earth, given as 


= 


Red - fg) 


v/l - (2f_ - fh cos^ 6 


(3-87) 


and r = magnitude of the satellite position vector 

R^ = equatorial radius of the Earth 

fg = Earth's flattening coefficient 

6 = declination of the satellite (it is assumed that 6 
equals the geocentric latitude of the subsatellite 
point) 

If the density is assumed to be maximum at the apex of the 
bulge, then the cosine variation between the maximum and 
minimum density profiles is 


P- (h) = p . (h) 

' min ' 


* (2) (3-88) 
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where y is the angle between the satellite position vector 
and the apex of the diurnal bulge. 

The cosine function in Equation (3-88) can be determined 
directly as 


cos 


n 


1 = 

2 


1 + cos Y 


] • [i 


u 


1 n/2 


B 


2r 


(3-89) 


where = satellite position vector (TOE coordinates) 

= unit vector directed toward the apex of the 
diurnal bulge (TOE coordinates) 

For FEDS, n has been assigned the value of 6. 

Vector Uq has the following components; 


X 

= cos 


cos (a 

s 

+ X) 


= cos 


sin (a 

s 

+ X) 

“b 

= sin 

«s 




(3-90) 


where 6g = declination of the Sun 

Os = right ascension of the Sun 

X = lag angle between the Sun line and the apex of 
the diurnal bulge (approximately 30 degrees) 

3.4.4 SOLAR RADIATION PRESSURE ACCELERATION 

The model for acceleration ~a, due to direct solar radiation 
acting on a satellite, is 




m 



(3-91) 
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where v = an eclipse factor such that 

V = 0 if the satellite is. -in shadow {umbra) 

V = 1 if the satellite is in sunlight 

0 < V < 1 if the satellite is in penumbra 

The vector from the Sun to the satellite is given by 


r 


vs 



(3-92) 


where = mean solar flux at one astronomical unit, di- 

vided by the speed of light 

r“g = position vector of the Sun 

T" = position vector of the satellite 

^SUN “ astronomical unit 

C„ = solar radiation pressure constant 

I\ 

A = surface area of the satellite 

m = mass of the satellite 

r = magnitude of the vector 'T 
vs vs 

The solar radiation pressure constant is given by 

= 1 + n (3-93) 

where n is the reflectivity of the surface. All of the 
factors listed above, except , are constant for a given 
satellite . 

A simple cylindrical shadow model is used to determine the 
eclipse factor. From Figure 3-1, it is apparent that the 
satellite is in sunlight (v = 1) if the vector product 
('r‘ • U^) is greater than zero, where is the satellite 
position vector relative to the Earth, and Ug is the solar 
position unit vector relative to the Earth. 
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SUNLIGHT 



Figure 3-1. Cylindrical Shadow Model 

If this product is less than zero and the vector to the 
satellite along the normal to the solar unit vector 

= t' - (T • Ifg) iTg (3-94) 

has a magnitude less than the Earth's radius, then the satel- 
lite is in shadow (i.e., v equals 0); otherwise, it is as- 
sumed that the satellite is in sunlight and v equals 1. 

3.5 VARIATIONAL EQUATIONS 

The variational equations are differential equations describ- 
ing the rate of change of state parameters. The variational 
factors are used to solve the following equation (see Sec- 
tion 3.2.3) ; 

Y = A(t) Y + B(t) Y + C(t) (3-95) 
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8S21/B2 


where 


A(t) 


B(t) 

C(t) 





3r" 

at: 


ac 


D 


(3-96) 


To reduce computation in FEDS, the acceleration vector used 

in the variation equations is reduced to the following: the 

Earth's central body acceleration, the zonal terms 

X T , X , , X T in the nonspherical gravitational acceleration, 
^J2 J3 J4 ^ 

■a‘ 2 , and the atmospheric drag acceleration, *a£). Using this 

reduced acceleration model, A(t) can be expressed as 


A(t) 



3r‘ 3^^ 3"? 


(3-97) 


The Earth's central body acceleration is defined in Equa- 
tion (3-58) . The partial derivatives with respect to the 
position vector are given by 
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(3-98) 
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The nonspherical gravitation acceleration is defined in 
Equation (3-76) . The zonal terms in the nonspherical poten- 
tial, have m equal to zero, such that the partial de- 

rivatives of these terms with respect to the longitude, X, 
are zero. In this case. Equation (3-76) reduces to 



3r 

T 

m 

+ — ^ 
a<j> 

(frf 

(3-99) 

where the zonal 

potential , 

♦z' 

given 

by 



00 

= ^ Z 



P°(sin (})) 

(3-100) 


n=2 


and 

/ R \ 

(n + 1) C° P°(sin (j» 

n (3-101) 

-) =n = ♦> 

n=2 




n=2 


The partial derivatives of ^ with respect to r are ob- 
tained by differentiating Equation (3-99) as follows; 
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(3-103) 


The required partial derivatives of are as follows; 
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To minimize computations, the symmetry property of the sec- 
ond partial derivatives of utilized as shown below. 

These partial derivatives are obtained by differentiating 
Equation (3-100) as follows; 
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(3-104) 


n=2 
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(3-105) 
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The partial derivatives of r and <j> with respect to r are 
given in Equations (3-79) . The required second partial de- 
rivatives of r and (j) with respect to "r* are obtained by 
differentiating Equations (3-79) with respect to ~t, yield- 
ing 



(3-106) 



where dx/d'f', dy/d'r', and 9z/8‘r“ are (1, 0, 0)/ (0, 1, 0) , 
and (0, 0, 1), respectively. 

Taking into account the symmetry properties of the second 
partial derivatives of r, <(), and X 
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and the fact that the potential function, ip, satisfies 

2 

Laplace's equation, A \p = 0 ; therefore 



(3-108) 


reduces the amount of computation required. 

The atmospheric drag acceleration is defined in Equa- 
tions (3-83) and (3-84) . The partial derivatives of "aj^ 
are as follows: 
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The partial derivatives of the density with respect to 
tion are derived from Equation (3-88) as follows: 
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The partial derivative of the height with respect to ■? 
obtained by differentiating Equation (3-86) , yielding 
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where 
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SECTION 4 - DATA PROCESSING FOR TRANSPONDER INTERFACE 


This section contains algorithms used in processing data 
received from and to be sent to the transponder, including 
algorithms to convert the raw Doppler data message into 
engineering units and to create a frequency control word. 
Communication with the transponder is accomplished using an 
Intel 8086 microprocessor interface provided by Code 530. 
FEDS forms data to be sent to the transponder into a message 
and transmits it to the microprocessor via a standard RS232 
terminal port. The microprocessor then translates the mes- 
sage into pulses to be sent to the various connectors on the 
transponder. For data traveling from the transponder to 
FEDS, the microprocessor collects the pertinent information 
from the transponder ports, forms a message, and transmits 
to FEDS again via a standard terminal port. 

4.1 REDUCTION OF TDRSS DOPPLER DATA 

In FEDS, the only observation data collected are one-way 
TDRSS Doppler measurements. The raw measurement consists of 
a nondestruct Doppler count of a nominal bias frequency 
added to a Doppler sample over a fixed time interval. The 
count is cumulative because the accumulator is reset only 
between passes, not between successive measurements. 

The one-way Doppler measurement is performed by transmitting 
a signal from a ground transmitting station to a forward- 
link TDRS. The TDRS coherently translates the signal to the 
tracking frequency of the user transponder. The transponder 
then accumulates a Doppler measurement for transmission to 
FEDS. 

Although all measurements should be valid due to the reset- 
ting of the accumulator and interception by the microproc- 
essor interface of measurements with invalid lock flags, the 
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first measurement of a tracking pass is ignored. The rejec- 
tion of the first measurement should not significantly af- ■ 
feet the results of the estimation because there is expected 
to be an abundance of observation data. The following algo- 
rithm is used to preprocess all Doppler observations; 







17 

240 


(4-1) 


where f (t. ) = observed Doppler shift (hertz) averaged over 
^0 the time interval between t^_^^ and t^ 

= value of the accumulator (counts) at time t^ 

= change in the value of the ac- 
cumulator (counts) between t^.^ and t^ 

M = number of Doppler samples added to the ac- 
cumulator during the time interval between 
ti_i and ti (nominally 40000) 

fg * frequency bias (nominally 2^1 = 2097152) 

K = rational multiplier of the Doppler sample 
(nominally 1/4) 

fu - frequency of the oscillator used by the 
transponder in forming a Doppler sample 
(nominally 19.056392 megahertz) 

During a nominal pass, the Doppler accumulator is expected 
to overflow at least once. When overflow is detected, 
used in the preprocessing of the Doppler observa- 
tion, will be adjusted as follows; 


(4-2) 

4.2 CREATION OF FREQUENCY CONTROL WORDS 

The frequency control word is produced by FEDS and trans- 
mitted to the transponder to enable acquisition of the TDRSS 
signal. FEDS outputs an initial control word at a fixed 
timespan (nominally 60 seconds) prior to the beginning of a 
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tracking interval and will provide an updated control word 
at a fixed frequency (nominally one control word every 
6 seconds) until acquisition occurs or until the end of the 
scheduled tracking pass. Output of control words will re- 
sume if the tracking signal is lost before the scheduled end 
of the pass. 

The first bit of the frequency control word is not used by 
FEDS. The second bit is always set to 1. The remaining 
14 bits of the control word are a signed binary (two's com- 
plement) number computed as follows: 

f ,(t.) • 2^^ • 3121 

CW(t. ) (4-3) 

^ ^ref 

where CW(ti) = data portion of the frequency control word 

(-8196 £ CW £ 8195) for output at time ti 

^d(^i) ” predicted Doppler shift at time ti 

f^g£ = reference frequency 
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SECTION 5 - ESTIMATION LOGIC 


Reference 9 provides the basis for the following section, 
which describes the estimation technique used in FEDS: a 

sliding batch differential correction (SBDC) algorithm. The 
procedure uses a classical weighted least-squares method 
solved by Newton iteration. The SBDC moves in discrete 
steps along the observation data, performing a state estima- 
tion at each step (slide) . The output solve-for state for 
each slide becomes the a priori input state for the follow- 
ing slide. The major assumptions underlying this approach 
are as follows: 

1. Given moderate changes in the solve state from 
iteration to iteration, the measurement partial derivatives 
do not vary enough to affect the DC process. Thus, the par- 
tial derivatives are computed for the first iteration of 
each slide and are held constant for subsequent iterations 
at which moderate state changes occur. 

2. If the a priori state is reasonably close to the 
actual state, all of the spurious measurements can be edited 
during the first DC iteration. Thus, an iterative residual 
edit loop is always performed for the first DC iteration 
and, on subsequent iterations, is done only when the linear- 
ity constraint is violated. 

3. Idle time in the system should be used to advan- 
tage. The algorithm provides a precompute phase that ini- 
tiates the next DC slide before all of the measurement data 
is available. 

4. The solutions are not noise dominated; as a result, 
data can be sampled down to produce a reduced observation 
set that decreases the data storage requirements. 
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The following subsections describe the SBDC algorithm as 
implemented in FEDS; this algorithm resembles the estimator 
used in GTDS (Reference 2) . 

This description is not meant to be a rigorous derivation of 
the techniques used in FEDS. The first three subsections 
describe one complete slide of the estimator, and later sub- 
sections cover the preliminary edit criteria (Section 5.4) 
and the special-case logic used to complete a partially pre- 
computed slide (Section 5.5) . 

Each slide consists of one or more iterations through the 
measurement data to determine a state vector, X, that best 
fits the data in a least-squares sense. The state vector, 

X, that is solved for contains a maximum of 10 parameters 
(the first 6 are mandatory, whereas the last 4 are optional) 

X — [ r , r , Cjj, b ] 

where 'f' = user spacecraft position vector in inertial 
Cartesian coordinates 

r = user spacecraft velocity vector in inertial 
Cartesian coordinates 

Cq = coefficient of drag term 

b = user spacecraft clock error terms (second-order 
polynomial) 

Each iteration can be divided into three functional areas; 

• An initial summation of the normal matrix and asso- 
ciated batch estimation matrices 

• An inner iterative process that consists of the 
computation of the state correction vector followed 
by an adjustment to the matrices previously accumu- 
lated 

• End-of-iteration testing for slide convergence or 
divergence and linearity violations 
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Depending on the result, the slide may either terminate with 
a new solution (converge) or terminate with no new solution 
(diverge) or continue by performing another iterative cycle. 

5.1 INITIAL SUMMATION OF BATCH ESTIMATION MATRICES 

The first functional area of the iterative process performs 
one complete pass through the measurement data. The com- 
puted measurements from the observation model (Section 6) 
are used to form the measurement residuals, which are then 
subjected to a preliminary edit procedure (Section 5.4) . 
Partial derivatives are computed for the nonedited measure- 
ments, and the normal matrix is formed. To facilitate the 
description of this algorithm, the following standard index- 
ing has been adopted. 

Parameter Meaning 

Subscript i Iteration-dependent variable 

Subscript j Observation-record-dependent variable 

The initial summation is performed by processing each avail- 
able observation record in the following manner: 

1. Compute measurement residuals (o - c) : 

(o - c) j = Oj - Cj (5-1) 


where Oj = preprocessed observed measurement (Section 4) 
Cj = computed measurement (Section 6) 


2. Perform preliminary editing based on maximum allow- 
able residual test and, for new measurements only, check the 
measurement geometry (Section 5,4). 
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3. For the nonedited measurements, compute the partial 
derivatives, Sj, the observation at time tj, with respect to 
the solve-for variables, X, at epoch time t^; 



3f (t .) aT (t .) 

- o - JL - -^o 1 

3r(t.) 3r"(t.) 

3 3 




3C, 


3b 


(5-2) 




3X, 


3X, 




3X 

3X 


0 


0 


where = the observation equation used to compute c^ 

The first matrix on the right is explicitly determined from 
the observation equations in Section 6. The second matrix 
is obtained by integrating the variational equations (Sec- 
tion 3.5). The values of a^ are used to form a single row 
of the F-matrix. 


T 

4. The normal matrix, F WF, and two other estimation 
T 

matrices, F WF (o - c) and Q, are then updated with current 
measurement information as follows: 



a 


where a equals the measurement standard deviation. 
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f'^WF. . = F*^WF. . , + (aT * w * a.) 
i/J D D 


(5-4) 


f'^WF(o - c). . = F^WF(o - c). . , + a. * w * (o - c) . (5-5) 

^ r 3 ^ 1 3 3 3 


Q + (o - o; 

Q. . = ^ 

1/3 a 


(5-6) 


The total number of observations used in the summations, n, 

is also incremented. Only the upper half and diagonal of 

T 

the symmetric F WF matrix are stored. 

5.2 STATE CORRECTION AND INNER EDIT 

1. The correction to the state, AX., is formed by 

T — ^ *p 

solving the system (F WF)^ AX^ = F WF(o - c)^ as follows: 

Air. = (f\f)T^ F^WF(o - c) . (5-7) 

11 1 

T 

where F WF^ is inverted using the Schur identity method of 
partitioning (Reference 2, Section 8.6.1). 

2. Estimation statistics are computed for later use. 
The standard error fit is as follows; 


hi - 

aIT^ * f’^wf(o - C) ^ 

/ 

n - 1 


The residual root mean square is as follows: 


(5-8) 



(5-9) 
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3. Inner loop editing is performed on the first itera- 
tion or if the linearity constraint was violated on the pre- 
vious iteration (Section 5.3) as follows: 

Predicted residuals, rj, are computed and tested 

r. = ^ ((o - c). - (a.) (Air.)] (5-10) 

J g* j j 

Given m, an input scaling parameter, each rj is compared 

with m.S.. If the r . is greater, the observation is edited, 

T 

and Its contribution is subtracted from the F WF . , 

T ^ 

F WF (o - c)^, and matrices. The total number of obser- 
vations used, n, is also decremented. 

4. If no observations are edited, inner editing is 

terminated. Otherwise, new and RMS^ statistics are com- 
puted based upon the modified sums and a new state correc- 
tion vector, AX^, is obtained. Inner editing will terminate 
at this point if either the maximum number of inner loops 
has been performed or if < E, where E is the 

standard error fit input tolerance. Otherwise, the inner 
edit loop is repeated until one of the above conditions is 
satisfied. The last state correction vector computed during 
this process becomes the state correction vector, A)T^, for 
this iteration. 

5.3 END-OF-ITERATION TESTING 

Various tests are made at the end of each iteration to de- 
termine whether the solution state has been attained and, if 
not, whether the correction has violated linearity con- 
straints. The convergence/divergence tests are performed in 
the following order: 

1. Divergence indicators from previous operations are 
checked for one of three conditions: error return from in- 

version of the normal matrix, observation timespan less than 
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the minimum required, or all observations in the latest pass 
have been edited. 

2. Convergence checking is performed next. The slide 
converges if both of the following are true; 



(5-lla) 


(5-llb) 


or if 





RMS. - 


(5-llc) 




Si 

V L* V^IIV 2 

where conv^, 
erances. 

conv 2 / and conv^ are uplinked convergence 

tol- 

3. 

lowing 

The slide diverges if 
is true: 

i = 1 and either of the 

fol- 




1 AP^I 

> diVj^ 

(5-12a) 




1 AV. 1 
1 

> aivj 

(5-12b) 


where div^ and div 2 are uplinked divergence tolerances. 

The slide diverges if i > 1 and either of the following is 
true: 


lAP^I > m2 


(5-12C) 
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(5-12d) 


lAV^I > m2 

where m 2 is a scaling factor. 

4. If the maximum number of iterations has been per- 
formed, a test for reduced convergence is made by multiply- 
ing conVj^, conv 2 , and conv^ by 3 and checking for 
convergence based on the new tolerances. Otherwise, the 
slide diverges. 

5. If neither convergence or divergence has occurred, 
a linearity check is made and another complete iteration is 
performed. The linearity test determines whether or not new 
partial derivatives, a^, are to be computed and the inner 
loop edit performed during the next iteration. The 
linearity constraint is satisfied (i.e., the test succeeds) 
if lAr^l < e^^ and lAr^l < e 2 / where e^^ and 62 are uplink 
linearity tolerances. 

If the test just described succeeds, no new partial deriva- 
tives are computed, and the partial derivatives from the 

m 

previous iterations are used to compute the [F WF(o-c)] 
matrix while the normal matrix remains unchanged. 

5.4 PRELIMINARY EDIT CRITERIA 

During the initial processing of the measurement data de- 
scribed in Section 5.1, two edit checks are performed to 
detect and remove anomalous data from the estimation proc- 
ess. The first simply verifies that the measurement resid- 
uals, (o - c)j, fall within acceptable limits. The 
residual is compared with the maximum allowed, and the ob- 
servation is edited if the residual exceeds the maximum. 

The next edit check is based on the measurement geometry 
between the user spacecraft and the TDRS spacecraft. 
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For the geometry edit test, two parameters are computed; H, 
the height of the TDRS-to-user ray path segment, and C, the 
central angle formed between two vectors (the Earth center- 
to-TDRS spacecraft and the Earth center-to-user space- 
craft) . These parameters are compared with two uplink 
parameters, Hq and C^; the measurement is edited only if 
H < Hq and C > Cq, where H and C are defined by the fol- 
lowing ; 


H = IP I - R (3 > 90^) 

user e 


' "’usee' ® - ‘'e (6 < 90 ) 


(5-13) 


where 


6 = cos 


-1 


^^user * ^^user ~ ^TDRS^ ^ 

ll^ I 1^ - I 

user user TORS 


(5-14) 


where Puser “ user satellite position (appropriate geocentric 
coordinates) 

Ptdrs ~ appropriate TORS position (appropriate geo- 
centric coordinates) 

Rg = equatorial radius of the Earth 
NOTE ; 3 is in units of degrees with 0° < 3 < 180°. 


C 


cos 


-1 


rp • p 1 
^ TORS user^ 


[IP 


TORS 


P I ] 
user ^ 


(5-15) 


where 0° < C < 180° 

5.5 SLIDE PRECOMPUTATION 

During precomputation of the next slide, an arbitrary epoch 
is chosen by advancing the epoch of the previous slide by a 
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fixed At. Therefore, the measurement partial derivatives 
computed during the precompute phase of estimation are ref- 
erenced to that arbitrary epoch. When additional measure- 
ments are made available to the estimator and when the 
actual estimation process begins as described in Sec- 
tion 5.1, a new epoch may be chosen. If the epoch used 
during precomputation is earlier than the new observation 
timespan, the epoch for the slide is set to the latest ob- 
servation time. If this is done, the partial derivatives 
are no longer referenced to the proper epoch time. To ad- 
vance the partial derivatives to the new epoch, the measure- 
ment partial derivatives are postmultiplied by the state 
transition matrix, which maps the Cartesian state from the 
old to new epochs. The updated partial derivatives are then 
used during the estimation process. If the epoch used dur- 
ing precomputation is later than the final new observation, 
estimation is referenced to the precomputation epoch-. 
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SECTION 6 - OBSERVATION MODELS 


This section contains the mathematical specifications for 
the observation models in FEDS. Section 6.1 contains the 
observation model for one-way TDRSS Doppler measurements. 
Section 6.2 contains the algorithm for computing the partial 
derivatives of the measurements with respect to the target 
satellite state. Section 6.3 contains miscellaneous models 
required to solve the observation and observation partial 
derivative equations. 

6.1 MEASUREMENT EQUATIONS 

A monochromatic source of electromagnetic radiation in 
motion with respect to the observer experiences a phenomenon 
referred to as Doppler shift. It can be approximated by the 
following relation; 


f 


ref c + V 


rel 


( 6 - 1 ) 


where f = observed, Doppler-shifted frequency 

fj-gf = transmitted, reference frequency (no relative 
motion between source and observer) 

Vrel = line-of-sight component of the relative velocity 
between the source and the observer 

c = speed of light 

For FEDS, observation measurements have the geometry shown 
in Figure 6-1. The figure shows a series of nodes that are 
connected by legs. The nodes are the end points of the sig- 
nal path (i.e., TDRSS spacecraft. White Sands Ground Track- 
ing Station (WSGT) , FEDS at GSFC, and a fictitious target 
spacecraft) ; the legs are the signal paths between the 
nodes. The arrows on the legs represent the direction of 
signal propagation. 
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(TORS) 
NODE 2 



Figure 6-1. Signal Geometry for Doppler Observation 

During the demonstration, WSGT will transmit at a frequency 
given by 


f,(t - 


6t. 


- 6t^) = 


(t 


- «t3) 


02 (t 


- 6t, + 6t,) ref 


( 6 - 2 ) 


where f^. ( t) 
fref 

a. (t) 


frequency of signal transmitted at time t 

demonstration reference frequency (nominally 
2287.5 megahertz) 

Doppler shift due to leg i at receive time t 
time of signal propagation for leg i 
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It should be noted that leg 1 and leg 3 signal propagation 
times are affected by transmission through the atmosphere. 

The propagation time is given by the following relation: 

6ti = + ARC^ (6-3) 


where Aili = time of signal propagation in a vacuum for 

leg i (light-time correction; see Section 6.3.1) 

ARCi = refraction correction associated with leg i, 
applied only for ground-to-space or space-to- 
ground legs (see Section 6.3.2) 

WSGT will compute Cj^(t) based on predicted fictitious space- 
craft ephemeris, actual ground antenna position, and real- 
time TORS solutions. Equation (6-3) indicates that the 
signal would arrive at a fictitious target spacecraft at a 
nearly constant frequency. 

The signal will arrive at the experiment transponder at a 
frequency given by 


fg(t) = aj^(t - a^{t) f^(t - 6t3 - 6t^) (6-4) 


where fQ(t) is the frequency received by the transponder 
located at GSFC. Combining Equations (6-2) and (6-4) , the 
received frequency becomes 


(t) 

^G^^^ ^ a2(t - 6t^ + 6t2) ^ref 


(6-5) 
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Substituting Equation (6-1) into Equation (6-5) twice and 
dropping negligible terms, the instantaneous Doppler shift 
is obtained as 


fo(t) " ^ref 


'ref 


(P2 - P 


( 6 - 6 ) 


where fQ(t) = the instantaneous Doppler shift at time t 


dp. 


Pi = 


i dt 


(6-7) 


The various leg distances, pj^, are defined as follows: 


P 3 = l^(t) - r^(t - St^) I 


( 6 - 8 a) 


P 2 = - 6tj + 6t2> - ■?‘^(t - St^) I (6-8b) 


where 'r‘^(t) = position vector of the transponder at time t 

^^(t) = position vector of TDRSS satellite at time t 

^f(t) = simulated position vector of the fictitious 
target satellite at time t 

Averaging over the interval At, the Doppler observation is 
computed as 


/' 


fp(t) dt 


Fj^(t) = 


t-At 


At 


ref 

cAt 


Ap_ - Ap- 


(6-9) 


where the change in leg length, Apj^, is defined as 


Ap^ = p^(t) - p^(t - At) 


( 6 - 10 ) 
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6.2 PARTIAL DERIVATIVES OF MEASUREMENTS WITH RESPECT TO 
TARGET STATE 


In general, the partial derivatives of an observation are 
computed as follows: 


3f (tj^) 
3lT3(to) 


3 X 3 (tR - 6t3 + 6t2) 


$3 (tj^ - 6t3 + 6t2, t^) 

( 6 - 11 ) 


where f(t) = a measurement at time tag t 

= state vector of node 3 at time t 
$ 3 (t, t^) = 3lT3 (t)/3lT3 (t^) = state transition matrix 
t^ = epoch time 
t^ = observation time tag 

The Doppler measurement partial derivatives are obtained 
from Equation (6-9) as follows: 


/ ^^^2 \ 

3?3(to)’=“ I 


( 6 - 12 ) 


Since leg 3 is independent of the target state. 


3Ap. 


3X3 (t^) 


= 0 


Equation (6-12) reduces to 


3Ap, 


8X3 (t^) \ax3(t^) 


(6-13) 


(6-14) 


6-5 


9681 



where 


3Ap. 


3X3(t^) 


3P2(V 


3Xj(tj^ - + 6t2) 

0 

0 

0 


3P2(tj^ - At) 


* f(t - 6t3 + 6t2, t^) 


(6-15) 


3X^(tj^ - At - 6t3 + 6t2) 


0 

0 

0 


* *(tj^ - At - 6t3 + 6t2, t^) 


with 


3P2(t) 


r^(t - 8t^ + 6t^) - r-^(t - 6t^) 


ax^ct - + 6t2) 


(6-16) 


6.3 FREQUENCY BIAS EFFECTS ON OBSERVATION MODELING 

In FEDS, frequency biases are modeled as a second-order 
polynominal given by 


AF(t) = AF(tg) + AF(tg)(t - t3) + AF(tg)(t - t^) ^ (6-17) 


where AF(t) = frequency bias at time t 
t = epoch time 

AF(t^) = frequency drift at epoch 

• • 

AF(tg) = frequency acceleration at epoch 

When frequency biases are included in the transmitted 
frequency. Equation (6-5) becomes 


fQ ( t ) 


Q2 (t 


Oto ( t) 

h + 6t2) ^ref ot3(t)aj^(t - 6t3) 
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Since a^Ct) and aj^(t - 6t^) are near unity, the Doppler 
shift can be obtained from Equation (6-6) as follows; 


f_(t) = (p, - pO + AF(t) 


(6-19) 


The averaged Doppler over the interval At is then 


FD(t) = (AP 2 - AP 3 ) + AF(t) - (6-20) 


Partial derivatives with respect to frequency biases can be 
obtained by differentiating Equation (6-20), as follows; 


A / 1_ \ 

3Fo(t) 

3AF(t^) 

3AF(tg) 

6.4 MISCELLANEOUS MODELS REQUIRED FOR OBSERVATION MODELING 

The miscellaneous models required to compute both the 
observation and the observation partial derivatives are 

• Newton-Raphson light-time corrector 

• Tropospheric refraction corrections 


(t - t^)^ - (t - At - tg)^ 
6 At 


(6-23) 


(t - tg)^ - (t - At - tg)2 
__ 


( 6 - 22 ) 


6-7 


9681 



6.4.1 NEWTON- RAPHSON LIGHT-TIME CORRECTOR 


The Newton-Raphson iterative method for calculating the time 
at a transmitting node, given the receiving node's time, 
(backward light-time correction) is 




c 


+ 




1 - 


)-nl 

hn) 

c ■ 

- ' (^T„) 



(6-24) 


where 


t' = new estimate of the transmission time from 
■^n node n 



■R 


n+1 


r (t) 
n ' 


n , n+1 



(t) 


previous estimate of the transmission time 
from node n /initialized to time t_ \ 

\ 

receive time at node n+1 
position vector at time t 

unit vector between nodes n and n+1 (leg n) 
velocity vector of node n at time t 


c = speed of light 


Ignoring negligible terms. Equation (6-24) can be reduced to 



(6-25) 


or 


6t 


N 






c 


(6-26) 


6-8 


9681 



p - p ' 
N N 


The equation is solved N times until 
6.4.2 TROPOSPHERIC REFRACTION CORRECTIONS 


< meter. 


The general approach to the problem of refraction correction 
of measurements is to obtain the most accurate corrections 
possible that are consistent with rapid computations. To 
preserve the raw tracking data in its original form, refrac- 
tion corrections are applied to the computed values of the 
measurements. Furthermore, the corrections are applied only 
to range. The Doppler observation is corrected by relating 
it to corrected range values. The method of computing tro- 
pospheric refraction presented here was taken from Refer- 
ence 10. 


Measurement Corrections . The refraction-corrected, computed 
measurements are obtained as follows: 


R 


cc 


R 


C 


+ 


(6-27) 


where R_^ = corrected, computed range 

CC 

Rc = computed distance between two nodes (uncorrected 
range) 

AR^ = refraction correction 

The Doppler shift is corrected for refraction effects by 
using the corrected range in Equation (6-9) . 

Tropospheric Correction . Tropospheric refraction correc- 
tions to range are assumed to be functions of two varia- 
tions: elevation angle and surface ref ractivity . Range 

correction is tabulated for nominal surface ref ractivity , 
indexed by elevation angle.- The tabulated values were com- 
puted by a ray tracing algorithm presented in Reference 10. 

The standard value of surface refractivity was chosen to be 
340 N units, because it is the average value of the range of 
surface ref ractivities encountered at the stations (280 N to 
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400 N units) . The approximation for range correction is in 
terms of variation of the surface refractivity from the 
standard. The method, then, is an algorithm that produces 
tropospheric refraction correction to range, and these re- 
fraction corrections agree with ray tracing results to 
within the specified design criteria. 

The refraction correction to range is therefore obtained by 
evaluating 



(6-28) 


where AR«ji = interpolated range correction for the 

standard surface refractivity 

3ARfp/3N = first-order partial derivatives of the 
range correction with respect to change 
in surface refractivity from the standard 

3^ARii>/3n 2 = second-order partial derivative of the 
range correction with respect to change 
in surface refractivity from the stand- 
ard; the superscript bar denotes values 
for the standard refractivity obtained 
through interpolation 


The method of computing the tropospheric refraction correc- 
tions for the range proceeds as follows. For the specified 
elevation angle, the necessary parameters are obtained from 
the tabulated data by numerical spline interpolation. (The 
spline interpolation method is the mathematical analog of 
the draftsman's mechanical spline, which is a long, very 
flexible, slender device used to pass a smooth curve through 
many data points. The technique is presented in Refer- 
ence 11.) Equation (6-28) is evaluated to obtain the re- 
fraction correction to the range. 
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